May27
Devon DeRaad
5/26/2020
convert p values to qvalues
#read in all phase1 significance values:
all.sigs<-read.csv(file="~/Desktop/anoph.3.march.2020/all.sigs.csv")
table(all.sigs$chrom)##
## 2L 2R 3L 3R X
## 638059 946059 524268 737401 198852#convert 0 bayescenv q-vals to 1e-4, to avoid infinite values when plotting with -log10()
all.sigs$prec.q[all.sigs$prec.q == 0]<-.0001
all.sigs$temp.q[all.sigs$temp.q == 0]<-.0001
all.sigs$hum.q[all.sigs$hum.q == 0]<-.0001
#convert lfmm prec p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
qvals<-qvalue(all.sigs$prec.p[all.sigs$chrom== i])
all.sigs$prec.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
#convert lfmm hum p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
qvals<-qvalue(all.sigs$hum.p[all.sigs$chrom== i])
all.sigs$hum.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
#convert lfmm temp p values to genome-wide FDR adjusted q values
for (i in levels(all.sigs$chrom)){
qvals<-qvalue(all.sigs$temp.p[all.sigs$chrom== i])
all.sigs$temp.p[all.sigs$chrom== i]<-qvals[["qvalues"]]
}
rm(qvals) #space saver
#give order for plotting manhattan plots
all.sigs$order<-c(all.sigs$pos[all.sigs$chrom== "2R"],
all.sigs$pos[all.sigs$chrom== "2L"]+61545105,
all.sigs$pos[all.sigs$chrom== "3R"]+61545105+49364325,
all.sigs$pos[all.sigs$chrom== "3L"]+61545105+49364325+53200684,
all.sigs$pos[all.sigs$chrom== "X"]+61545105+49364325+53200684+41963435)lfmm manhattans
#LFMM precip plot -log10(signficance vals) as manhattan plots
plot1<-ggplot(data=all.sigs, aes(x=order, y=-log10(prec.p), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.prec.outliers<- all.sigs[all.sigs$prec.p < .01,]
nrow(lfmm.prec.outliers) #check number of outliers## [1] 17113#LFMM hum plot -log10(signficance vals) as manhattan plots
plot2<-ggplot(data=all.sigs, aes(x=order, y=-log10(hum.p), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.hum.outliers<- all.sigs[all.sigs$hum.p < .01,]
nrow(lfmm.hum.outliers) #check number of outliers## [1] 10981#LFMM temp plot -log10(signficance vals) as manhattan plots
plot3<-ggplot(data=all.sigs, aes(x=order, y=-log10(temp.p), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
lfmm.temp.outliers<- all.sigs[all.sigs$temp.p < .01,]
nrow(lfmm.temp.outliers) #check number of outliers## [1] 12396#
grid.arrange(plot1,plot2,plot3, nrow=3)
#bayescenv manhattans
#bayescenv precip plot -log10(signficance vals) as manhattan plots
plot1<-ggplot(data=all.sigs, aes(x=order, y=-log10(prec.q), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant precip outliers
baye.prec.outliers<- all.sigs[all.sigs$prec.q < .01,]
nrow(baye.prec.outliers) #check number of outliers## [1] 3857#bayescenv hum plot -log10(signficance vals) as manhattan plots
plot2<-ggplot(data=all.sigs, aes(x=order, y=-log10(hum.q), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant hum outliers
baye.hum.outliers<- all.sigs[all.sigs$hum.q < .01,]
nrow(baye.hum.outliers) #check number of outliers## [1] 3688#bayescenv temp plot -log10(signficance vals) as manhattan plots
plot3<-ggplot(data=all.sigs, aes(x=order, y=-log10(temp.q), col = chrom))+
geom_scattermore(pointsize = 1.2)+
theme(panel.grid=element_blank(), panel.background=element_blank())+
geom_hline(yintercept=-log10(.01), linetype="dashed")
#pull all significant temp outliers
baye.temp.outliers<- all.sigs[all.sigs$temp.q < .01,]
nrow(baye.temp.outliers) #check number of outliers## [1] 10195grid.arrange(plot1,plot2,plot3, nrow=3)
par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(na.omit(lfmm.hum.corrs$id[abs(lfmm.hum.corrs$diff) < .05]))
sig.pos<-sig.pos[1:27]
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(na.omit(lfmm.hum.corrs$id[abs(lfmm.hum.corrs$diff) > .8]))
sig.pos<-sig.pos[1:27]
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%check out the correlations of allele frequency patterns with humidity
dev.off()## null device
## 1#visualize
plot(lfmm.hum.corrs$diff)
abline(h=0, col="red")
plot(lfmm.hum.corrs$phase1, lfmm.hum.corrs$phase2)
hist(lfmm.hum.corrs$diff)
plot(-log10(lfmm.hum.outliers$hum.p), all.cases$diff)
plot(-log10(lfmm.hum.outliers$hum.q), all.cases$diff)
plot(lfmm.hum.outliers$hum.fst, all.cases$diff)different classes of relationships?
dev.off()## null device
## 1layout(matrix(c(1,1,1,1,2,3,4,5),2, 4, byrow=TRUE))
plot(lfmm.hum.corrs$phase1)
abline(h=c(.4,.1,-.25), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 > .4])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < .4 & lfmm.hum.corrs$phase1 > .1])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < .1 & lfmm.hum.corrs$phase1 > -.25])
abline(v=c(-.2,.2), col="red")
hist(lfmm.hum.corrs$diff[lfmm.hum.corrs$phase1 < -.25])
abline(v=c(-.2,.2), col="red")strongly positive correlation
par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 > .4])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%weakly positive correlation
par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < .4 & lfmm.hum.corrs$phase1 > .1])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%weakly negative correlation
par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < .1 & lfmm.hum.corrs$phase1 > -.25])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%strongly negative correlation
par(mfrow=c(3,3))
#plot
sig.pos<-as.vector(lfmm.hum.corrs$id[lfmm.hum.corrs$phase1 < -.25])
sig.pos<-sample(sig.pos, size=27, replace = FALSE)
j<-1 #initialize loop tracker
pb <- txtProgressBar(min = 0, max = length(sig.pos), style = 3) #initialize progress bar##
|
| | 0%for (i in sig.pos){
frq<-as.numeric(as.vector(phase1.all.freqs[phase1.all.freqs$id == i,3:16]))
frq2<-as.numeric(as.vector(phase2.all.freqs[phase2.all.freqs$id == i,3:23]))
plot(sort.pops.phase1$hannual,frq, main = i, ylim =c(0,1))
abline(lm(frq~sort.pops.phase1$hannual))
points(sort.pops.phase2$hannual,frq2, col="red")
abline(lm(frq2~sort.pops.phase2$hannual), col="red")
setTxtProgressBar(pb, j)
j<-j+1
}##
|
|=== | 4%##
|
|===== | 7%##
|
|======== | 11%##
|
|========== | 15%##
|
|============= | 19%##
|
|================ | 22%##
|
|================== | 26%##
|
|===================== | 30%
##
|
|======================= | 33%##
|
|========================== | 37%##
|
|============================= | 41%##
|
|=============================== | 44%##
|
|================================== | 48%##
|
|==================================== | 52%##
|
|======================================= | 56%##
|
|========================================= | 59%##
|
|============================================ | 63%
##
|
|=============================================== | 67%##
|
|================================================= | 70%##
|
|==================================================== | 74%##
|
|====================================================== | 78%##
|
|========================================================= | 81%##
|
|============================================================ | 85%##
|
|============================================================== | 89%##
|
|================================================================= | 93%##
|
|=================================================================== | 96%
##
|
|======================================================================| 100%